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Abstract 

We present a novel approach to investigating the divergence history of biomes 
and their component species using single-locus data prior to investing in multilo- 
cus data. We use coalescent-based hierarchical approximate Bayesian computation 
(HABC) methods (MsBayes) to estimate the number and timing of discrete diver- 
gences across a putative barrier and to assign species to their appropriate period of 
co-divergence. We then apply a coalescent-based full Bayesian model of divergence 
(IMa) to suites of species shown to have simultaneously diverged. The full Bayesian 
model results in reduced credibility intervals around divergence times and allows 
other parameters associated with divergence to be summarized across species as- 
semblages. We apply this approach to 10 bird species that are wholly or patchily 
discontinuous in semi-arid habitats between Australia's southwest (SW) and south- 
east (SE) mesic zones. There was substantial support for up to three discrete periods 
of divergence. HABC indicates that two species wholly restricted to more mesic 
habitats diverged earliest, between 594,382 and 3,417,699 years ago, three species 
from semi-arid habitats diverged between 0 and 1,508,049 years ago, and four di- 
verged more recently, between 0 and 396,843 years ago. Eight species were assigned 
to three periods of co-divergence with confidence. For full Bayesian analyses, we 
accounted for uncertainty in the two remaining species by analyzing all possible 
suites of species. Estimates of divergence times from full Bayesian divergence models 
ranged between 429,105 and 2,006,355; 67,172 and 663,837; and 24,607 and 171,085 
for the earliest, middle, and most recent periods of co-divergence, respectively. This 
single-locus approach uses the power of multitaxa coalescent analyses as an efficient 
means of generating a foundation for further, targeted research using multilocus 
and genomic tools applied to an understudied biome. 



Introduction 

Environmental change during the Plio-Pleistocene had a ma- 
jor influence on the composition and distribution of species 
globally (e.g., Hewitt 1996, 2004; Avise and Walker 1998; 
Provan and Bennett 2008). Northern hemisphere ice sheets 
made vast regions uninhabitable and caused the majority of 
species there to persist only in southern refugia (Provan and 
Bennett 2008 ). In contrast, Australia was not glaciated during 
this period, but its biota still faced environmental transforma- 
tions associated with glacial cycles. The continent's Miocene 
mesic climates and environments (e.g., Kershaw et al. 1994; 



Martin 2006) became the more xeric ones that dominate the 
continent today (McLaren and Wallace 2010). Pleistocene 
glacial cycles caused repeated expansion and contraction of 
the Australian arid zone culminating in the most severe arid 
phase at the Last Glacial Maximum (20-18 kA) (Williams 
et al. 2009). Consequently, some species and lineages once 
presumably more widespread across southern Australia are 
now restricted to the southwest (SW) and southeast (SW) or 
east of the continent. A sample of such histories is our focus 
in this paper. 

Relatively few Australian biogeographic regions and bar- 
riers have been studied using comparative phylogeography 
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(Avise 2000). Key examples, however, mostly apply to eastern 
Australian tropical and sclerophyllous forest environments 
(e.g., Joseph et al. 1995; Hugall et al. 2002; Schauble and 
Moritz 2001; Stuart-Fox et al. 2001; Garrick et al. 2008). 
Comparative phylogeographies of arid and semi-arid zone 
bird species have revealed variable levels of phylogeographic 
structure, some species lacking structure (Joseph and Wilke 
2007) and others showing clear signals of historical isola- 
tion (Joseph and Wilke 2006; Joseph and Omland 2009; 
Kearns et al. 2009). Concerning studies of major phylogeo- 
graphic barriers to gene flow in particular, the Carpentarian 
Barrier in northern Australia (Jennings and Edwards 2005; 
Lee and Edwards 2008) and the Black Mountain Corridor 
in the Wet Tropics rainforest (Dolman and Moritz 2006) 
have additionally been studied with multilocus approaches. 
Together, phylogeographic studies of single species, barri- 
ers, and regions, whether from single locus or multilocus 
data, are bringing new dimensions to understanding Aus- 
tralia's Plio-Pleistocene biogeography (Byrne 2008; Byrne et 
al. 2008, 2011). Nonetheless, some biogeographic regions 
have been under-represented. In these regions, a single-locus 
comparative phylogeographic approach using the power of 
multitaxa coalescent analyses is an efficient way of provid- 
ing sound foundational knowledge for future multilocus 
studies. 

Limitations of mitochondrial DNA (mtDNA) as an infer- 
ential tool in phylogeography are well known. They arise from 
stochastic variation in the coalescent process (Kuo and Avise 
2005), selection (Ballard and Whitlock 2004), and sex bias 
in fitness or dispersal behavior (Hare 2001). These limita- 
tions can be overcome by surveying more loci (Edwards and 
Beerli 2000) or more species using comparative phylogeog- 
raphy (Avise 2000). A remaining problem with comparative 
phylogeography is that signatures of disparate timing of di- 
vergence may reflect either continuous divergence, multiple 
discrete periods of divergence, or stochastic variation in the 
coalescent or differences in population demography around 
a simultaneous divergence (e.g., Dolman and Moritz 2006; 
Hickerson et al. 2006). Further, a single locus offers low preci- 
sion for estimating parameters such as population divergence 
time (Edwards and Beerli 2000; Hudson and Turelli 2003). 
A solution to this conundrum comes from recently devel- 
oped coalescent-based methods (Hickerson et al. 2006) that 
make it possible to re-examine the utility of mtDNA alone 
in testing for temporal congruence among species across a 
putatively common barrier. These methods can account for 
demographic variation and coalescent stochasticity and es- 
timate the number and timing of periods of co-divergence 
among a suite of species affected by the same barrier. Fur- 
ther, they can exploit variation between species with identical 
histories of divergence to estimate divergence time (akin to 
variation between loci within a single species divergence). 



Coalescent-based hierarchical approximate Bayesian com- 
putation (HABC) methods for estimating the number and 
timing of periods of co-divergence have shed light on a 
number of well-recognized barriers to gene flow (Beau- 
mont 2010). For simplicity, hereafter, these discrete episodes 
when multiple species co-diverged will be referred to as "co- 
divergence events." Two historical co-divergence events were 
inferred across both the Isthmus of Panama (Hickerson et 
al. 2006) and Baja California's mid-peninsula (Leache et al. 

2007) . Likewise, comparative mtDNA phylogeographies of 
hydrothermal species in the East Pacific Rise (Plouviez et al. 
2009), birds of the Queen Charlotte Islands (Topp and Winker 

2008) , and across the Isthmus of Tehuantepec (Barber and 
Klicka 2010) have also benefited from the statistical rigor 
of coalescent-based HABC using MsBayes (Hickerson et al. 
2006, 2007). Most recently, lives et al. (2010) demonstrated 
the advantages of HABC analyses in testing long-standing 
hypotheses of colonization or persistence of intertidal and 
subtidal invertebrate communities of the west and east coasts 
of the Atlantic Ocean. 

Here, we apply these methods to study the Plio-Pleistocene 
histories of 10 bird species that are either isolated in mesic 
zones of Australia's SE and SW or which are patchily discon- 
tinuous in semi-arid habitats between those regions. Several 
putative barriers to gene flow have been recognized in this 
broad biogeographic region (see Fig. 1, inset) (Schodde and 
Mason 1999). The Lake Eyre Basin/Eyrean Barrier (Keast 
1961; Ford 1974; Schodde 1982) has been affirmed by molec- 
ular data as a major phylogeographic break in Australian ring- 
neck parrots (Barnardius zonarius) (Joseph and Wilke 2006) 
and the splendid fairy-wren {Malurus splendens) (Kearns et 
al. 2009). Here, we focus attention on the Nullarbor Plain, 
which is another barrier of likely pivotal biogeographic im- 
portance because of its geographically central location in 
southern Australia (Fig. 1). The Nullarbor Plain is an exten- 
sive region of some 200,000 square kilometers (Webb and 
James 2006) of limestone karst, formed through elevation of 
limestone sediments about 14 million years ago (Mya) (Li et 
al. 2004) and continuously exposed since then (Sheard and 
Smith 1995). It is dominated by chenopod shrubs. Currently, 
it extends from the continent's south coast northwards to 
the southern fringes of the more heavily vegetated central 
Australian arid zone (see review in McKenzie and Robinson 
1987). Coupled with its fringing xeric woodlands and shrub- 
lands, it forms the core of a strong biogeographic barrier of 
approximately 700 km between Australia's temperate mesic 
zones in the SW and SE (Fig. 1). At glacial maxima, lower 
sea levels allowed rivers to drain more effectively, resulting 
in fewer swamps and floodplains (Haq et al. 1987; Macphail 
et al. 1994). This, together with regression of the shoreline 
causing subaerial exposure of the Great Australian Bight (Li 
et al. 2004), offered opportunity for some species to disperse 
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Figure 1. Map of Australia showing locations 
of Eyrean and Nullarbor Barriers in relation to 
other putative, minor barriers in southern 
Australia (after Schodde and Mason 1999), 
and sampling localities for 10 bird species. 
Mesic species: white circles = Petroica 
boodang, white squares = Phylidonyris 
novaehollandiae, white diamonds = 
Melithreptus lunatus; semi-arid zone species: 
gray diamond = Gliciphila melanops, black 
square (white outline) = Drymodes 
brunneopygia, gray square = Glossopsitta 
porphyrocephala, black x = Malurus 
pulcherrimus, black diamond (white outline) = 
Eopsaltria griseogularis, black cross = Malurus 
splendens, gray circle = Climacteris rufus. 
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and recolonize, potentially allowing gene flow between the 
SW and SE. 

Crisp and Cook (2007) used relaxed molecular clock dat- 
ing to estimate the time of interspecific divergences in sister 
species pairs of plants on either side of the Nullarbor Plain. 
Divergence times between 23 endemic SW plant lineages and 
sister lineages endemic to SE were mostly relatively ancient 
(clustering around 13-14 Mya). Mid-Miocene aridification 
and/or elevation of the Nullarbor Plain as a vicariant bar- 
rier were inferred as the driving agents of these divergences. 
Three younger divergences of 2-4 Mya were thought to be 
due to severe drying around this time (Fujioka et al. 2005). 
Among birds, phylogeographic studies of Australian mag- 
pies (Gymnorhina tibicen) (Toon et al. 2007), southern emu- 
wrens (Stipiturus malachurus) (Donnellan et al. 2009), and 
musk ducks (Biziura lobata) (Guay et al. 2010) have also 
revealed divergent monophyletic intraspecific lineages on ei- 
ther sides of the Nullarbor Plain and its fringing semi-arid 
habitats. 

In this study, we first test for temporal congruence of di- 
vergence across the Nullarbor Plain using HABC. For this, 
we use the MsBayes software pipeline that allows for varia- 
tion in stochasticity of the coalescent process and variation 
in demographic histories among taxa (Hickerson et al. 2006, 
2007). We estimate the number of co-divergence events that 
generated the current distributions of the 10 bird species, as- 
sign suites of species to the appropriate co-divergence event 
and estimate their timing (Hickerson et al. 2006, 2007). 



Our approach is summarized in Figure 2 and as follows. 
For the species assemblage divergence(s) uncovered using the 
HABC approach, we apply a full Bayesian approach (IMa) 
(Hey and Nielsen 2007) to estimate population demographic 
parameters associated with the divergence of suites of species 
(hereafter termed species assemblage divergence). To our 
knowledge, this is the first time IMa has been used in a species 
assemblage context. IMa (Hey and Nielsen 2004, 2007) is 
intended to fit an Isolation-with-Migration (IM) model to 
genetic data from multiple loci sampled from a single pair of 
divergent populations. IMa's goal is to estimate demographic 
parameters associated with the two populations' divergence 
from their ancestral population (Nielsen and Wakeley 2001; 
Hey and Nielsen 2004). Rather than estimating population- 
splitting parameters for the SW and SE population pair of 
each species using multiple loci, we use IMa for a novel 
purpose. Having first used MsBayes to identify discrete co- 
divergence events relevant to 10 southern Australian bird 
species, we then used IMa to estimate demographic param- 
eters associated with each co-diverging species assemblage. 
Thus, for a given co-divergence event, rather than our ap- 
proach and data exhibiting stochasticity of the coalescent 
process among loci, they represent stochasticity across mul- 
tiple species affected by the same co-divergence event. IMa 
allows inheritance scalars (h) to be estimated as parameters in 
the model, such that h for each locus is free to vary during the 
course of the Markov chain as a function of the data and the 
model. There are two biological reasons for allowing h to vary 
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Figure 2. A flow chart outlining the sequence of HABC and Full-Bayesian analyses used in estimating parameters associated with the divergence of 
species assemblages. Box colors are indicative of the type of analysis: HABC (MsBayes) = blue; full Bayesian (IMa) = yellow; comparison and inference 
= green or red). 



among loci in IM (Hey and Nielsen 2004): first to account 
for biases in the effective sex ratio of males and females, and 
second the effective population size of particular loci could 
be reduced relative to other loci in the genome through re- 
current selective sweeps (Maynard Smith and Haigh 1974; 
Gillespie 2000) or background selection (Charlesworth et al. 
1993). In the context of running IMa as a taxon assemblage 
divergence model, invoking the option to estimate h allows 
the model to account for differences (within reason) in ef- 
fective population sizes among 10 SW and SE population 
pairs of southern Australian birds. In this context, estimating 
h allows effective population sizes to differ between species 
rather than between separate loci within a species. A limita- 
tion, however, is that effective population sizes of ancestral 
and the two divergent populations are not able to vary id- 
iosyncratically with respect to each other across all species. 
Effective population sizes of each of SW and SE population 
pair and their ancestral population are thus assumed to vary 
consistently across all taxa. This may be a considerable as- 



sumption if species have responded to habitat fluctuations 
differently. 

Materials and Methods 
Sampling 

We sampled four southern Australian mesic zone bird 
species that have populations isolated in SW and SE (tax- 
onomy and nomenclature follows Christidis and Boles 
2008): scarlet robin (Petroica boodang), new holland hon- 
eyeater (Phylidonyris novaehollandiae), white-naped hon- 
eyeater (Melithreptus lunatus), and tawny-crowned hon- 
eyeater (Gliciphila melanops). We sampled six species of 
semi-arid southern Australia, each having populations iso- 
lated wholly or partially in SW and SE: purple-crowned 
lorikeet ( Glossopsitta porphyrocephala), southern scrub-robin 
(Drymodes brunneopygia), blue-breasted fairy-wren (M. pul- 
cherrimus), western yellow robin (Eopsaltria griseogularis), 
splendid fairy-wren (M. splendens), and rufous treecreeper 
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(Climacteris rufus), the last two being patchily continuous 
in semi-arid and arid habitat to the north of the Nullarbor 
Plain. Two hundred seventy-eight individuals (n = 10-48 
per species) were sampled from these 10 species. This in- 
cluded 15 newly acquired specimens of M. splendens from 
the Great Victoria Desert north of the Nullarbor Plain, which 
had been identified in that species as a major sampling gap, 
and 32 ND2 sequences from M. splendens (Kearns et al. 2009) 
downloaded from GenBank (EU144272-EU144303). To fo- 
cus attention on the Nullarbor Plain and its fringing semi-arid 
habitats, we excluded samples east of the Eyrean Barrier (Fig. 
1) if Fst for a given species was significant across that barrier. 
Museum specimen and locality details of all individuals stud- 
ied, and mtDNA genealogies are provided in supplementary 
material. 

DNA extraction and sequencing 

Total cellular DNA was extracted according to the salting-out 
method of Miller et al. (1988). We amplified the ND2 
gene region of the mtDNA genome using primers L5216 
(5'-GGCCCATACCCCGRAAATG-3') and H6313 (5'- 
ACTCTTRTTTAAGGCTTTGAAGGC-3') (http://people. 
bu.edu/msoren/primers.html). Polymerase chain reactions 
(PCR) were performed in 25 [A reactions containing 1 x 
reaction buffer (Eppendorf), 0.2 mM dNTP, 0.2 fiM each 
primer, 1 Unit HotMaster™ Taq DNA Polymerase (Eppen- 
dorf), which specifically has no requirement for additional 
MgCh. Cycling conditions included an initial denaturation 
of 94°C for 2 min followed by a touch-down protocol of 34 
total cycles of 94°C denaturation for 30 sec, 62°C-58 °C 
annealing for 30 sec, 65°C for 45 sec, and a final extension of 
65° C for 5 min. PCR products were purified with Ampure 
magnetic beads (Agencourt), and sequenced directly in both 
directions on an ABI 3700 PRISM DNA analyzer. 

Summary statistics and tests of assumptions 

DNA divergence indices (nucleotide diversity, n; nucleotide 
divergence, Da) were calculated using DnaSP v5.0 (Librado 
and Rozas 2009). Addressing concerns about selective neu- 
trality of avian mtDNA diversity (Zink 2005), we tested for 
neutrality using MK tests (McDonald and Kreitman 1 99 1 ) in 
DnaSP v5.0 (Librado and Rozas 2009). 

Analyses of divergence of species 
assemblages 

A flow chart (Fig. 2) outlines the sequence of HABC and 
Full-Bayesian analyses used in estimating parameters associ- 
ated with the divergence of species assemblages. Individual 
analyses are described in more detail below. 
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Tests for temporal congruence of divergence 

Hierarchical Bayesian model tests under an ABC framework 
implemented in the MsBayes software pipeline (Hickerson 
et al. 2007) tested for temporal congruence of divergence 
times among populations in SE and SW. The MsBayes soft- 
ware pipeline simultaneously estimates three hyper parame- 
ters that quantify the variability in divergence times of the 10 
SW and SE population pairs: ( 1 ) mean population divergence 
time (E[r ]) in units of fi per generation, (2) degree of discor- 
dance (£2) as the variance of population divergence time r, 
divided by the mean of r across all population pairs, such that 
the greater the variance of population divergence time rela- 
tive to the mean divergence time across all population pairs, 
the less support there is for simultaneous divergence, and (3) 
number of discrete co-divergence events (*) (Hickerson etal. 
2006). A vector of observed summary statistics is calculated 
from DNA sequence data of each SW and SE population pair. 
Random draws of the hyper-parameters and sub-parameters 
from a coalescent model of population divergence with or 
without postdivergence gene flow (see Hickerson et al. 2007 
for details) were used to simulate finite sites DNA sequence 
data from 10 population pairs. In the model, ancestral and 
descendent population sizes are free to vary, thus allowing for 
the general model of isolation by colonization. From these 
simulations a vector of summary statistics was produced. 
From 5,000,000 iterations, acceptance/rejection with regres- 
sion algorithms were used to compare the observed summary 
statistics vector with the summary statistics vector from the 
simulated data. Polychotomous regression was used to es- 
timate and local linear regression was used to estimate 
Q and E(t). Hyper-posteriors were constructed using 1800 
accepted draws. 

We used the default summary statistics (n, net it, Wat- 
terson's 0, and 1/Tajima's D) for acceptance/rejection. We 
tried a range of upper limits for sub-parameters 6 (current 
population diversity per site per generation) and ancestral 
population size (relative to current population size) and set- 
tled on those that best fit the data (9 upper = 0.025, ancestral 
population size upper limit = 0.25). We initially tried high 
upper limit of r = 10. Given the low r estimates from these 
trials, we reduced the upper bound of r upper to 3.0. Anal- 
yses were run with migration rate upper limit set to 0.0, but 
allowing some migration at migration rate upper limit = 0.5 
did not affect results. Recombination rate upper limit was set 
to 0.0 (assuming no recombination in mtDNA). 

HABC analyses were initially performed on the 10 SW 
and SE population pairs allowing * to vary. For the full 
Bayesian analysis that follows (see Flowchart; Fig. 2), it is 
important that the number of discrete co-divergence events 
is not underestimated, or else multiple divergences will be 
analyzed as a single event. Overestimation of the number of 
discrete co-divergence events is less important, except in that 
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it may reduce the power to estimate divergence times. For this 
reason, the Bayes Factor (BF) was calculated to test the level 
of support for 4* < "highest posterior estimate" versus * > 
"highest posterior estimate." BFs were also calculated for the 
model of simultaneous divergence (*4> = 1 vs. * > 1) and 
a model in which diversification is continuous rather than a 
discrete number of divergence events (i.e., via rare dispersal 
events across pre-existing geographic or climatic barriers; ^ 
<10 vs. * = 10). For tests of simultaneous divergence £2 
with local linear regression has been shown to be a reliable 
indicator (Hickerson et al. 2006). Therefore, BFs were also 
calculated for £2 = 0 versus £2 > 0.0. The level of support 
was determined using the scale suggested by Jeffreys (1961). 
HABC analyses were then performed with * fixed to the 
highest posterior estimate of * from the initial analysis. To 
test the confidence in our assignment of species to each co- 
divergence event, we ran subsets of taxa allowing * to vary. If 
the hyper-parameter estimate of £2 is very close to zero and * 
= 1, there is support for a single co-divergence event for that 
subset of taxa. BFs were calculated for ^ = 1 versus *)/ > 1 to 
give an indication of the level of support (Jeffreys 1961). 

Estimating demographic parameters 
of community divergence 

IMa (Hey and Nielsen 2004, 2007) and our approach to its use 
for divergence of species assemblages are described above (see 
Introduction). To reiterate briefly, the goal is to estimate de- 
mographic parameters associated with the divergence of two 
daughter populations from an ancestral population (Hey and 
Nielsen 2004; Nielsen and Wakeley 2001). The parameters 
are scaled by fi, the geometric mean across loci of muta- 
tion rate/year/locus and are effective population diversity for 
ancestral population, 9 A, and pairs of daughter populations 
since divergence 9 1 and 92(9 = 4 N/x), directional migration 
rates, ml and m2, and population divergence time, t. 

The data were divided into datasets representing discrete 
species assemblage divergences as uncovered with MsBayes, 
and each species assemblage divergence was analyzed sep- 
arately. Infinite sites and HKY are the only models of nu- 
cleotide evolution available for DNA sequence data in the 
IMa program. We used the HKY model because mtDNA loci 
are assumed to be nonrecombining and the data do not fit 
the infinite sites model. Preliminary metropolis Monte Carlo 
Markov chain simulations (20 million steps) were required to 
optimize the number of chains, heating model, heating pa- 
rameters, and upper bounds for the parameter priors. Broad 
parameter limits were subsequently refined. Six chains with 
geometric heating (gl = 0.8, g2 = 0.92) were found to main- 
tain a sufficient level of swapping between chains (between 
chains 4 and 5 was 5-79% and between chains 0 and 1 was 
81-92%). Two final simulations for each divergence were 
carried out with different random seeds, 50 million steps, 
and burn-in of five million. Chain swapping and genealogy, 
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branching, and mutation rate mixing were sufficient, there 
were no trends present in parameter trend plots and effec- 
tive sample sizes were ample (minimum ESS per run were 
between 72 and 523). Two final runs per analysis confirmed 
convergence upon comparable parameter distributions. 

Converting divergence time estimates to a 
demographic scale 

To convert divergence time estimates to a demographic scale, 
we applied an imported mutation rate based on a general 
avian divergence rate of 2.1% (Weir and Schluter 2008). The 
general avian mutation rate was generated with a time re- 
versible (GTR+r) model and was cross-validated with 74 
avian molecular clock calibrations generated from fossil and 
biogeographic events across 12 taxonomic avian orders and 
over the last 12 Ma (Weir and Schluter 2008). Divergence 
time estimates from MsBayes were calculated using popula- 
tion divergence time in years, t = (E[r ]) x (6Ave//x), where 
(E[t ]) is the mean population divergence time in units of /x 
per generation, 9 Ave is half of the upper bound of the 9 prior 
used in the MsBayes analysis, and /x is the per site per gener- 
ation mutation rate. We assumed an average generation time 
of two years based on published data (Higgins 1999; Higgins 
et al. 2001; Higgins and Peter 2002) and acknowledge that 
divergence time estimates will be sensitive to assumptions of 
generation time. Divergence time estimates from IMa were 
calculated using population divergence time t = t //x, where t 
is the divergence time parameter estimate from IMa, and \± is 
per locus per year mutation rate. These divergence times were 
compared to those calculated using the imported mutation 
rate applied directly to average net divergence (Da) estimates. 
General patterns of effective population sizes for SW, SE, and 
ancestral populations are described for each co-divergence 
event. 

Results 

Two hundred forty-six ND2 sequences were collected from 
the 10 species. ND2 sequences did not contain stop codons or 
indels, consistent with not having amplified pseudogenes. All 
sequences were submitted to GenBank (JQ27369-JQ27614; 
see Table SI for individual level identification). 

Summary statistics and tests of assumptions 

Haplotype diversity ranged from 0.739 to 0.963 (Table 1). 
Nucleotide diversity (n) within SW or SE populations varied 
by an order of magnitude, from 0.049% in E. griseogularis 
(SW) to 0.54% in G. melanops (SW) (Table 2). Five of the 
10 population pairs are reciprocally monophyletic, having 
Da between 0.19% and 4.35% (Table 2). The other five have 
net divergences up to three orders of magnitude less, be- 
tween 0.001% and 0.08%. Three of these show no structure 
between SW and SE populations and individuals from the 
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Table 1. Polymorphism and haplotype statistics based on ND2 sequences. 





Sequence 


No. of 


No. of 


No. 


of parsimony No. of 


Haplotype diversity 


Species 


length 


individuals 


polymorphic sites informative sites haplotypes 


HD (+SD) 


Melithreptus lunatus 


1033 


30 


63 


55 


20 


0.963 














(0.945-0.981) 


Petroica boodang 


1020 


38 


45 


32 


1 8 


0.913 














(0.887-0.939) 


Phylidonyris novaehollandiae 


1033 


29 


33 


20 


1 6 


0.936 














(0.909-0.963) 


Malurus splendens 


987 


48 


38 


1 6 


22 


0.75 














(0.682-0.818) 


Climacteris rufus 


1022 


32 


23 


9 


1 7 


0.925 














(0.897-0.953) 


Gliciphila melanops 


865 


27 


28 


16 


18 


0.943 














(0.911-0.975) 


Eopsaltria griseogularis 


I uzz 


z / 


7 


4 


5 


U. / oJ 














(0.732-0.834) 


Drymodes brunneopygia 


865 


14 


10 


3 


9 


0.934 














(0.889-0.979) 


Malurus pulcherrimus 


1022 


10 


9 


1 


6 


0.778 














(0.641-0.91 5) 


Glossopsitta porphyrocephala 


1024 


23 


13 


2 


12 


0.739 














(0.638-0.84) 


Standard deviations (SD) are given in parentheses. 










Table 2. Nucleotide diversity within and nucleotide divergence between SW and SE populations based on ND2 sequences. 




No. of 


No. of 












individuals 


individuals 


Reciprocally 


7T West 


jt East 


e-WWest (9 -W East Da (%) 


Species west 


east 


monophyletic? 


(%) (SD) 


(%) (SD) 


(SD) (SD) 


(SD) 


Melithreptus 1 4 


16 


Yes 


0.22 


0.40 


0.0027 0.0050 


4.35 


lunatus 






(0.18-0.26) 


(0.34-0.46) 


(0.0014-0.0041) (0.0029- 


-0.0071) (3.73-4.96) 


Petroica boodang 1 0 


28 


Yes 


0.13 


0.18 


0.0017 0.0038 


2.85 








(0.09-0.17) 


(0.16-0.21) 


(0.0007-0.0027) (0.0023- 


-0.0053) (2.35-3.35) 


Phylidonyris no- 1 0 


19 


Yes 


0.27 


0.33 


0.0038 0.0033 


1.28 


vaehollandiae 






(0.21-0.34) 


(0.30-0.37) 


(0.0019-0.0056) (0.0019- 


-0.0048) (1.04-1.53) 


Malurus 1 4 


34 


Yes 


0.47 


0.14 


0.0052 0.0012 


0.98 


splendens 






(0.42-0.51) 


(0.08-0.20) 


(0.0039-0.0095) (0.0033- 


-0.0071) (0.72-1.24) 


Climacteris rufus 6 


26 


Yes 


0.34 


0.25 


0.0039 0.0036 


0.19 








(0.28-0.40) 


(0.22-0.28) 


(0.0017-0.0060) (0.0021- 


-0.0050) (0.10-0.29) 


Gliciphila 1 1 


16 


No 


0.54 


0.52 


0.0060 0.0085 


0.08 


melanops 






(0.44-0.64) 


(0.45-0.59) 


(0.0033-0.0087) (0.0051- 


-0.0119) (0.00-0.20) 


Eopsaltria 1 2 


15 


No 


0.05 


0.0013 


0.0010 0.0012 


0.04 


griseogularis 






(0.03-0.07) 


(0.11-0.15) 


(0.0003-0.0016) (0.0005- 


-0.0019) (0.00-0.07) 


Drymodes 3 


1 1 


No 


0.16 


0.27 


0.0016 0.0032 


0.03 


brunneopygia 






(0.08-0.23) 


(0.24-0.31) 


(0.0003-0.0029) (0.0016- 


-0.0048) (0.00-0.14) 


Malurus 2 


8 


No 


0.12 


0.0021 


0.0012 0.0030 


0.00 


pulcherrimus 






(0.06-0.18) 


(0.13-0.30) 


(0.0000-0.0024) (0.0014- 


-0.0046) (0.00-0.12) 


Glossopsitta por- 9 


14 


No 


0.23 


0.06 


0.0036 0.0012 


0.00 


phyrocephala 






(0.19-0.28) 


(0.04-0.08) 


(0.0018-0.0054) (0.0005- 


-0.0020) (0.00-0.05) 


Nucleotide diversity is given as n 


(as a percentage) and Watterson's 8 (0-W). Divergence is given as net divergence, Da (as a percentage). Standard 



deviations (SD) are given in parentheses. 
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Table 3. Tests of population expansion. 



Species 


Ramos and Onsin's R 2 West 


Ramos and Onsin's R 2 1 


iast Fu's F s west 


Fu's F s east 


Melithreptus lunatus 


0.113 


0.104 


-1.756 


-6.906 




P = 0.093 


P = 0.106 


P = 0.133 


P = 0.002 


Petroica boodang 


0.143 


3.855 


-1.547 


-7.819 




P = 0.102 


P = 0.001 


P = 0.091 


P = 0.001 


Phylidonyris novaehollandiae 


0.120 


0.128 


-3.796 


-0.578 




P = 0.055 


P = 0.425 


P = 0.014 


P = 0.922 


Malurus splendens 


0.072 


0.083 


-5.995 


-6.139 




P = 0.001 


P = 0.1 16 


P = 0.004 


P= 0.015 


Gimacteris rufus 


0.099 


0.093 


-2.941 


-4.186 




P = 0.0504 


P = 0.159 


P = 0.053 


P = 0.025 


Gliciphila melanops 


0.130 


0.076 


0.523 


-4.793 




P = 0.177 


P = 0.004 


P = 0.591 


P= 0.019 


Eopsaltria griseogularis 


n 1 aa 


U. I JO 


? 1 1A 
— Z . I Z 


— U.OU D 




P = 0.147 


p = 0.394 


P = 0.025 


p= 0.222 


Drymodes brunneopygia 


0.471 


0.121 


NA 


-2.403 




P = 0.310 


P = 0.064 




P = 0.054 


Malurus pulcherrimus 


0.500 


0.173 


NA 


-0.965 




P= 1.000 


P = 0.331 




P = 0.216 


Glossopsitta porphyrocephala 


0.084 


0.113 


-5.197 


-3.143 




P = 0.003 


P = 0.003 


P = 0.004 


P = 0.003 


Ramos-Onsins and Rozas's Rj 


and Fu's Fs are given with P values (a = 


; 0.05) based on coalescent simulations (1 0,000 replicates). P values that reached 


statistical significance are in bold. NA, not applicable owing to small 


sample size. 






Table 4. Observed summary statistics from 10 co-distributed SW and SE population pairs used for acceptance/rejection in MsBayes analyses. 


Species 


IT 


Netjr 


Watterson's 9 


1/Tajima's D 


Melithreptus lunatus 


0.024787 


0.041916 


0.015394 


0.036908 


Petroica boodang 


0.012806 


0.027713 


0.0105 


0.032003 


Phylidonyris novaehollandiae 


0.009034 


0.012496 


0.008135 


0.026547 


Malurus splendens 


0.007556 


0.010777 


0.008051 


0.027008 


Climacteris rufus 


0.002571 


0.002367 


0.003542 


0.014701 


Gliciphila melanops 


0.005678 


0.000836 


0.008516 


0.026703 


Glossopsitta porphyrocephala 


0.001258 


0.000396 


0.00344 


0.016245 


Eopsaltria griseogularis 


0.001143 


0.000298 


0.001777 


0.012167 


Malurus pulcherrimus 


0.002312 


0 


0.00376 


0.012275 


Drymodes brunneopygia 


0.002602 


0 


0.003686 


0.014251 



SW and SE share haplotypes (D. brunneopygia, G. porphy- 
rocephala, M. pulcherrimus) (Table 2). The remaining two 
species (E. griseogularis, G. melanops) have similarly unstruc- 
tured mtDNA diversity but there are no shared haplotypes 
among SW and SE individuals. Four SW populations and 
six SE populations showed signals of population expansion 
that reached statistical significance according to Fu's Fs and 
Ramos-Onsin and Rozas' R2 statistics (Table 3). MK tests did 
not reject neutrality. 

Tests for temporal congruence of divergence 

Observed summary statistics from the 10 SW and SE popu- 
lation pairs used in MsBayes analyses are presented in Table 
4. Estimates of Q, and * indicated that a model of three co- 



divergence events was best supported by the data (£2 = 0.394; 
Fig. 3). The BF for the model of simultaneous divergence (* 
= 1 vs. * > 1) is 2.059, or "barely worth a mention" (Jeffreys 
1961). Using threshold values of £2 = 0 versus £2 > 0.0, the 
BF = 0.273, indicating negative support. The BF for <P < 3 
versus *)7 > 3 is 8.189, indicating a substantial level of support 
for no more than three co-divergence events. Notably, for * 
< 5 versus * > 5, the BF is 135.585, indicating a decisive 
level of support for no more than five divergence events. The 
continuous model of divergence is overwhelmingly rejected; 
* = 10 versus * < 10, BF = 3.9 x 10" 10 . 

Fixing * = 3 (the highest posterior estimate) allowed 
MsBayes to estimate the timing of each species assemblage 
divergence and to estimate the number of species that di- 
verged at each co-divergence event. Results indicated that four 
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Figure 3. Posterior distributions of temporal 
concordance. (A) Prior Distribution (dashes) 
and posterior densities (solid line) of the 
number of co-divergence events, *. (B) 
Three-dimensional diagram of the joint 
posterior distribution of the degree of 
discordance, Q (the variance of r, divided by 
the mean of t across all population pairs), and 
the mean r across all population pairs (in units 
of neutral mutation rate, n) from MsBayes 
analyses. 




SW and SE population pairs diverged more recently (mode 
r = 8.297E-06), three diverged at mode x = 0.232, and two 
diverged relatively early (mode r = 1.145) (Table 5). The 
number of taxa assigned to each co-divergence event sums to 
nine (not 10). It is unclear whether the uncertainty lies at the 
boundary between the first and second co-divergence events, 
or between the second and third divergences. Accordingly, 
two of the 1 0 species ( fifth and eighth when ranked by net ji ) 
cannot confidently be assigned to their corresponding species 
assemblage divergence. Gliciphila melanops could have di- 
verged in either the most recent or the middle co-divergence 
event. Phylidonyris novaehollandiae could have diverged in 
either the deepest or the middle co-divergence event. Sub- 
sets of species were formed to account for this uncertainty in 
species assemblages for each co-divergence event (Table 5). 
Initially, E. griseogularis, D. brunneopygia, M. pulcherrimus, 
and G. porphyrocephala are confidently associated with the 
most recent co-divergence event. Gliciphila melanops could 
be associated with this event or the middle co-divergence 



event. Similarly, M. splendens and C. rufus are confidently 
assigned to the next oldest co-divergence event but P. novae- 
hollandiae could be associated with this event or the oldest 
co-divergence event. Finally, M. lunatus and P. boodang are 
confidently assigned to the oldest event. We accounted for the 
uncertainty in the placement of G. melanops and P. novaehol- 
landiae by analyzing all possible placements separately (see 
Table 5). To obtain an indication of the level of support for 
each possible suite of species for each co-divergence event, 
MsBayes analyses were repeated on all possible subsets with 
* unconstrained to estimate Q, and * (Table 5). Estimates of 
£2 very close to zero and * = 1 indicate positive support, and 
BFs were calculated to indicate the robustness of the support 
(Table 5). Parameter estimates of Q, were very close to zero 
for all but one subset of taxa where the estimate of £2 (0.091) 
was a significant outlier (P < 0.05) (Table 5). The highest 
posterior probability for * was equal to 1 in all but one ad- 
ditional subset of taxa (Table 5). This was when M. splendens 
and C. rufus were analyzed alone. BFs for VP = 1 versus ^ 
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Table 5. Number and timing of divergence events estimated from MsBayes analyses, assignment of species to each divergence event accounting for 


uncertainty, and level of support for a single divergence event for each possible suite of species. 






Co-divergence Number of Divergence Time (r) Divergence time (in years) 


Possible suites 


a 1> 


Event taxon pairs Mode 2.5%-97.5% Mode 2.5% -97.5% 


of species 


var(r)/E(r) 


lau i w — 4 o.zy/b-ub u.uuu— u.ddd iu u— Dyo,o4D 


Eopsdltrid 


n nnnn i 
U.UUUU 1 




griseogularis, 






Drymodes 






brunneopygia, 






Malurus 






pulcherrimus, 






Glossopsitta 






porphyrocephala 




(or b) 


4- Gliciphila melanops 


n nnnri 1 
U.UUUU I 


Tau 2 * = 3 0.232 0.000-1.267 276,430 0-1,508,049 


2 Malurus splendens, 


0.0130 3 2 




Climacteris rufus 






4- Gliciphila melanops 


3 n no in 1 

u.uyiu 1 




4- Phylidonyris 


U.UU4I I 




novaehollandiae 




1 (or 4) 


+ Gliciphila melanops, 


0.0180 1 




Phylidonyris 






novaehollandiae 




Tau 3 4i = 2 1.145 0.499-2.871 1,363,215 594,382-3,417,699 


Melithreptus lunatus, 


0.0006 1 




Petroica boodang 






+ Phylidonyris 




1 (or3) 


novaehollandiae 


0.0020 1 





1 A total of eight taxa were assigned to respective divergence events with confidence. Accounting for uncertainty, the number of taxon pairs for each 
divergence event is in parentheses. 

2 Malurus splendens and Climacteris rufus were assigned with confidence to the second divergence event, and were analyzed alone, however the 
divergence event includes either Gliciphila melanops or Phylidonyris novaehollandiae, or both of these species (Psi = 3 (or 4)). 
3 Q very close to zero and * = 1 indicates support for a single divergence event. A single divergence event for Malurus splendens and Climacteris 
rufus alone was not supported (* = 2). A single divergence event for Malurus splendens, Climacteris rufus and Gliciphila melanops was not well 
supported, with Q = 0.091 a significant outlier (P < 0.05). 



> 1 indicate that the assignment of species to the most recent 
divergence (with or without G. melanops), is the only co- 
divergence event indicating strong support for simultaneous 
divergence of its corresponding suites of species (BF = 27.8 1 7, 
29.598, respectively). BFs for the remaining suites of species 
indicate support is weak. For this reason, other combinations 
of species were analyzed (see Flowchart, Fig. 2), for example 
not including both M. splendens and C. rufus. There was no 
increased support for these additional species combinations, 
possibly reflecting reduced power of analyzing only two or 
three species in a single species assemblage divergence. 

Estimating demographic parameters of 
taxon assemblage divergence 

IMa analyses were performed on 10 SW and SE population 
pairs, divided into these same subsets of taxa. Divergence 
time estimates from IMa analyses are compared with diver- 
gence time estimates from MsBayes analyses and average net 
divergence in Table 6. 



IMa peak posteriors and 95% credibility intervals were 
plotted alongside the MsBayes modes and 95 percentile (Fig. 
4). Results were broadly consistent, but IMa results had sub- 
stantially reduced credibility intervals. Peak migration rate 
estimates were effectively zero in all but one scenario of co- 
divergence. Some migration from SW to SE was inferred 
when G. melanops was assigned to the first (most recent) co- 
divergence event (peak posterior ml = 0.825; 95% credibility 
interval: 0.122-3.662). 

Overarching patterns of effective population sizes of an- 
cestral, SW, and SE populations differed among the three co- 
divergence events. For the most recent co-divergence event, 
SW and SE assemblage population sizes were similar if G. 
melanops was included in the species assemblage, but if G. 
melanops was omitted the population size for the SW as- 
semblage was very large and could not be estimated. In 
this instance, the model of divergence is in doubt and the 
corresponding divergence time should not be relied upon 
without the comparative context provided here. Estimates 
of ancestral population size were two orders of magnitude 
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Table 6. Divergence time estimates from I Ma analyses for all suites of species are compared with those calculated using net divergence (Da) estimates. 





Possible suites 


Support for 


Divergence Time 










of species for species 


suite of 


V) 


Divergence time (in years 


) 






assemblage 


species 


IMa Analyses 


IMa Analyses 






Co-divergence 


divergence 


(MSBAYES 


Peak posterior 


Peak Posterior 


MsBayes Analyses 


Calculated from 


event 


analyses 


analyses) 


(2.5%-97.5%) 


(2.5%-97.5%) 


Mode (2.5%-97.5%) 


average Da 


Tau 1 


Eopsaltria 


Strong 


0.58 


56,954 


10 


9286 




griseogularis, 




(0.276-0.964) 


(25,607-89,439) 


(0-396,843) 






Drymodes 














brunneopygia, 














Malurus 














pulcherrimus, 














Glossopsitta 














porphyrocephala 












ir = 4 or 5 


+ Gliciphila melanops 


Strong 


0.956 


94,881 




1 5,048 








(0.284-1.844) 


(26,349-171,085) 






Tau 2 


1 Malurus splendens, 


Negative 


1.164 


111,119 


276,430 


279,048 




Climacteris rufus 




(0.724-4.060) 


(67,172-376,685) 


(0-1 ,508,049) 




i/c — 3 or 4 


\-\-\JIILipi Hid 


Negative 


u. / 00 


70 qqo 

/ 0| JJO 




I 30, / J5U 




melanops ) 














-\-Phylidonyhs 


VVedK 


1 


119. Q£3 




OQQ COO 




novaehollandiae 




(1.604-5.748) 


(148,818-533,296) 








+Gliciphila melanops, 


Weak 


1.615 


157,471 




301,786 




Phylidonyris 




(1.035-7.155) 


(96,027-663,837) 








novaehollandiae 












Tau 3 


Melithreptus lunatus, 


Weak 


9.375 


869,807 


1,363,215 


1,713,333 




Petroica boodang 




(4.625-21.625) 


(429,105-2,006,355) 


(594,382-3,417,699) 




f= 2 or 3 


+Phylidonyris 


Weak 


9.875 


914,267 




1,345,873 




novaehollandiae 




(4.875-16.125) 


(452,300-1,496,068) 







Figure 4. Comparison of divergence time 
estimates from HABC (MsBayes) and 
Full-Bayesian (IMa) analyses: MsBayes mode 
and 95 percentile and IMa peak posterior and 
95% Credibility intervals for each 
co-divergence event. Note that Tau 2, IMa (1), 
and IMa (2) are not shown as they had 
negative support for * = 1 and Q not close to 
0, respectively. 
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Analysis 



lower than current population sizes (for both SW and SE 
with G. melanops and only SE without G. melanops). For 
the middle co-divergence event, the ancestral population 
sizes were lower than the current population sizes, but by 
a more moderate scaling (2.7-5.7 times lower). Current SW 



and SE assemblage population sizes are similar. For the ear- 
liest co-divergence event, ancestral population size is com- 
parable to the current population size of the SE assem- 
blage, but the current SW population size is 1.8-4 times 
lower. 
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Discussion 

We compared demographic histories of 10 southern Aus- 
tralian bird species with discontinuous or patchily continu- 
ous distributions across southern Australia. Our broad aim 
was to explore how Plio-Pleistocene history of environmental 
change in southern Australia shaped evolution of the region's 
biota. We also wished to develop a foundation for targeting 
multilocus-based approaches to the same issue. Our specific 
aim was to use a coalescent-based HABC method, (MsBayes: 
Hickerson et al. 2006, 2007) to test for simultaneous diver- 
gence across one of the putative barriers in temperate south- 
ern Australia, in this case the Nullarbor Plain and its fringing 
semi-arid habitats (Fig. 1). Our central finding is that the 
model that is best supported by our data is one in which 
there were three separate co-divergence events in the history 
of the 10 SW and SE population pairs. The level of support 
for this model was substantial. There was only weak sup- 
port for a model of simultaneous divergence and a model 
of continuous divergence was overwhelmingly rejected. As 
the distributions of our focal species are either wholly dis- 
junct or patchily continuous across arid habitats north of the 
Nullarbor, we recognize the need to emphasize the distinction 
between co-divergence events and which barriers to gene flow 
are associated with them. The three temporal co-divergences 
may have been associated not just with the Nullarbor Plain 
but also with barriers other than the Nullarbor Plain and its 
fringing semi-arid habitats (e.g., Eyrean Barrier; see Intro- 
duction). Issues arising from the simple conclusion are now 
addressed. 

Regarding the assignment of species to particular co- 
divergence events, it is of interest to note that the two assigned 
to the earliest divergence event (M. lunatus, P. boodang) are 
restricted to more mesic habitats (predominantly wetter open 
forest and woodland formations). These are the habitats that 
would have been isolated longest across southern Australia 
(Byrne 2008; Byrne et al. 201 1). Indeed, our results are con- 
sistent with and offer some independent affirmation of Toon 
et al.'s (2010) conclusion that the eastern and western sub- 
species of M. lunatus that were recognized when this work 
began actually are two non-sister taxa, M. lunatus in the SE 
and M. chloropsis in the SW. Phylidonyris novaehollandiae 
may have diverged early with P. boodang and Melithreptus 
species, or in the middle divergence event. This species cer- 
tainly inhabits open forest and woodland, but is also common 
in drier heaths. The remaining seven species are more tolerant 
of drier, semi-arid habitats. 

Broad consistency of divergence time estimates derived 
from MsBayes and IMa is promising. This is all the more so 
because MsBayes is specifically designed for testing simulta- 
neous divergence of multiple taxa across a common barrier, 
whereas IMa is designed for multilocus inference within a sin- 
gle population split. To our knowledge, this is the first time 
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IMa has been applied to an assemblage of species. What made 
it possible was that MsBayes can first identify the different 
temporal co-divergence events prior to performing IMa on 
appropriate species assemblages. Assuming adequate perfor- 
mance of the inheritance scalar in accounting for different 
effective population sizes between SW and SE population 
pairs, inference from IMa analyses resulted in reduced cred- 
ibility intervals around divergence time estimates for each 
co-divergence event. As noted earlier, this approach assumes 
that effective population sizes of the ancestral and the two di- 
vergent daughter populations vary consistently with respect 
to each other across all species. The extent to which this as- 
sumption will affect inference from IMa analyses is yet to 
be tested. Caution is especially needed when drawing con- 
clusions on effective population sizes. In MsBayes analyses, 
these parameters are intentionally decoupled across taxa for 
this purpose. So the broad congruence of the divergence time 
parameter between MsBayes and IMa provides a level of 
confidence that this assumption is not greatly affecting the 
model of divergence of the species assemblage in the IMa 
analyses in this case presented here. 

Converting divergence time estimates to a 
demographic scale 

Converting divergence time to parameters on a demographic 
scale warrants some comment. The issue is contentious given 
uncertain estimations of mutation rates (Ho 2007). In addi- 
tion, applying a mutation rate generated from a phylogeny 
estimated using a GTR +r nucleotide substitution model 
to a population divergence model using a HKY nucleotide 
substitution model (as in IMa) is likely to underestimate di- 
vergence time (Arbogast et al. 2002). Given these limitations, 
the divergence time estimates have been converted to a de- 
mographic scale to make divergence time parameters broadly 
comparable across analyses and to provide only a very coarse 
timescale for these species assemblage divergences. 

According to the IMa analysis, which requires no assump- 
tion of generation time, the deepest divergence of two or 
three mesic SW and SE population pairs is from the middle 
to the early Pleistocene (95% CI: 429,105-2,006,355 years). 
Consistent with the potential for the HKY nucleotide substi- 
tution model to underestimate divergence time, especially in 
older divergences, this estimate is on the lower limit of Toon 
et al.'s (2010) divergence time for SW M. chloropsis and SE 
M. lunatus (2-4.6 Mya) from relaxed molecular-clock dat- 
ing. This divergence may be consistent with the progression 
of aridification following the step at 1.4-1.5 Mya identified 
by McLaren and Wallace (2010). In contrast, the most re- 
cent divergence of semi-arid SW and SE population pairs 
(95% CI: 25,607-171,085) is bound at the upper limit at the 
mid-Pleistocene, and includes the Last Glacial Maximum. 
Pertinent to this time period is that plant species with more 
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palatable leaves and fleshy fruits were present in the region of 
the Nullarbor between 180,000 and 400,000 years but are now 
restricted to remnant stands on the Nullarbor Plain's fringes 
(Prideaux et al. 2007). The limit for the most recent diver- 
gence is consistent with the floristic transition to present- 
day's dominant treeless chenopod shrubland habitats. Simi- 
larly, Hocknull et al. (2007) noted extinction of mesic fauna 
in central eastern Australia within this timeframe. 

Interspecific divergence times of 23 sister species pairs of 
plants endemic to either side of the Nullarbor Plain clustering 
around 13-14 Mya (Crisp and Cook 2007) are relatively an- 
cient compared to the three divergence times identified here 
in an assemblage of birds. Mid-Miocene aridification and/or 
elevation of the Nullarbor Plain as a vicariant barrier were in- 
ferred as the driving agents of the more ancient divergences. 
In Crisp and Cook's (2007) study, three younger divergences 
of 2-4 Mya are consistent with the deepest of the three di- 
vergence events in our study. This cluster of divergences was 
thought to be consistent with severe drying around this time 
(Fujioka et al. 2005). 

Synthesizing observations of population size from IMa re- 
sults, we suggest they are consistent with species being most 
dramatically affected prior to (or around the time of) the 
most recent co-divergence event (bound at the upper limit 
at the mid-Pleistocene, and including the Last Glacial Max- 
imum). This is consistent with the floristic transition to the 
present-day's dominant treeless chenopod shrubland habitats 
(discussed above). Perhaps not unexpectedly, the methodol- 
ogy used here supports this simple conclusion but in so doing 
suggests the utility of the method more generally. 

Despite the potential for gene flow in times of lower sea 
level during glacial maxima across land now submerged un- 
der the present day Great Australian Bight (Haq et al. 1987; 
Macphail et al. 1994; Li et al. 2004), there was generally no 
evidence for migration. Migration was inferred only for the 
most recent co-divergence event and indeed only when G. 
melanops was included. In this case, there was postdivergence 
gene flow from populations of the SW taxon assemblage to 
populations of the SE taxon assemblage. 

Conclusions 

Our central finding is that 10 bird species with discontinuous 
or patchily continuous distributions across southern Aus- 
tralia were differentially affected by three Plio-Pleistocene di- 
vergence events. The youngest two of the three co-divergence 
events are clearly within the Pleistocene, but the oldest po- 
tentially reaches back to the Pliocene according to the 95 
percentile of MsBayes results. It may be argued that we have 
not tested for the specific locations of these co-divergence 
events but the Nullarbor Plain and its fringing habitats are 
plausibly invoked as a barrier wholly or in part for each species 
assemblage divergence. The Eyrean Barrier to the east may 



also have contributed, especially to the earliest co-divergence 
event. This simple conclusion provides the basis for future 
studies on the biogeographic influence of Plio-Pleistocene 
environmental change in southern Australia, as well as for 
future multilocus-based approaches to the same issue. Our 
study highlights the synergistic benefits of jointly applying 
tests for simultaneous divergence with coalescent-based anal- 
yses of divergence. HABC methods (MsBayes) were first used 
to identify species assemblages shown to have simultaneously 
diverged, and then we applied a full Bayesian coalescent-based 
model of divergence (using IMa) in order to estimate parame- 
ters associated with co-divergence of this species assemblage. 
Divergence time estimates from MsBayes were broadly con- 
sistent with those from IMa applied to the divergence of 
species assemblages. Application of IMa in this context also 
enables exploration of migration and effective population 
sizes among co-diverging species assemblages. The joint ap- 
plication we have suggested of these methods utilizes the 
power of species assemblages to improve the insight that can 
be gained from single-locus data. 
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